/*
 * DenseLU.cs
 * 
 * Copyright (c) 2003-2005, dnAnalytics. All rights reserved.
*/

#region Using Directives
using System;
using dnAnalytics.Math;
using dnAnalytics.LinearAlgebra;
using dnAnalytics.LinearAlgebra.Native;
using dnAnalytics.Exceptions;
using dnAnalytics.Resources;
#endregion Using Directives

namespace dnAnalytics.LinearAlgebra.Decomposition
{
    /// <summary>
    /// The Compute method locks the object and call InternalCompute. InternalCompute will 
    /// decompose the matrix, overwriting the current matrix with the decomposition, and sets the
    /// the value of l1Norm, infNorm, and pivots.
    /// </summary>
    internal sealed class DenseLU : Algorithm, ILU
    {
        #region Fields

        private bool singular;
        private int[] pivots;
        private double det = Double.MinValue;
        private double l1Condition = Double.MinValue;
        private double infCondition = Double.MinValue;
        private double l1Norm;
        private double infNorm;
        private DenseMatrix matrix;

        #endregion Fields

        #region Constructors

        /// <summary>
        /// Constructs a DenseLU form an <see cref="Matrix"/> object.
        /// </summary>
        /// <param name="matrix">The matrix to use.</param>
        public DenseLU(Matrix matrix)
        {
            if (matrix == null)
            {
                throw new ArgumentNullException("matrix", Strings.NullParameterException);
            }

            if (matrix.Rows != matrix.Columns)
            {
                throw new NotSquareMatrixException(Strings.MustBeSquare);
            }

            this.matrix = new DenseMatrix(matrix);
        }

        #endregion Constructors

        #region Public Methods

        /// <summary>
        /// Returns a value indicating whether the matrix is singular.
        /// </summary>
        /// <returns><b>true</b> if the matrix is singular; otherwise, <b>false</b>.</returns>
        public bool IsSingular()
        {
            Compute();
            return singular;
        }

        Matrix ILU.LowerFactor()
        {
            return LowerFactor();
        }

        /// <summary>
        /// Returns the lower factor of the factorization.
        /// </summary>
        /// <returns>The lower factor of the factorization.</returns>
        public Matrix LowerFactor()
        {
            Compute();
            Matrix result = MatrixBuilder.CreateMatrix(matrix.Rows, MatrixType.Dense);
            LowerFactor(result);
            return result;
        }

        /// <summary>
        /// Places the lower factor in the result <see cref="Matrix"/>.
        /// </summary>
        /// <param name="result">The <see cref="Matrix"/> to store the lower factor in.</param>
        public void LowerFactor(Matrix result)
        {
            Compute();
            matrix.LowerTriangle(result);
            for (int i = 0; i < result.Rows; i++)
            {
                result.ValueAt(i, i, 1.0);
            }
        }

        Matrix ILU.UpperFactor()
        {
            return UpperFactor();
        }

        /// <summary>
        /// Returns the upper factor of the factorization.
        /// </summary>
        /// <returns>The upper factor of the factorization.</returns>
        public Matrix UpperFactor()
        {
            Compute();
            return matrix.UpperTriangle();
        }

        /// <summary>
        /// Places the upper factor in the result <see cref="Matrix"/>.
        /// </summary>
        /// <param name="result">The <see cref="Matrix"/> to store the lower factor in.</param>
        public void UpperFactor(Matrix result)
        {
            Compute();
            matrix.UpperTriangle(result);
        }

        /// <summary>
        /// The determinant of the matrix.
        /// </summary>
        /// <returns>The determinant of the matrix.</returns>
        public double Determinant()
        {
            Compute();
            if (singular)
            {
                return 0;
            }
            if (det == Double.MinValue)
            {
                det = 1.0;
                for (int j = 0; j < matrix.Rows; j++)
                {
                    if (pivots[j] != j)
                    {
                        det = -det * matrix.data[j * matrix.Rows + j];
                    }
                    else
                    {
                        det *= matrix.data[j * matrix.Rows + j];
                    }
                }
            }
            return det;
        }

        /// <summary>
        /// The L1 norm condition number of the matrix. The number is calculated as ||A|| * ||AInverse|| where ||.|| is the L1 norm.
        /// </summary>
        public double L1NormConditionNumber()
        {
            Compute();
            if (singular)
            {
                return Double.PositiveInfinity;
            }
            if (l1Condition == Double.MinValue)
            {
                l1Condition = matrix.lapack.ConditionNumber('1', matrix.Rows, l1Norm, matrix.data, pivots);
            }

            return l1Condition;
        }

        /// <summary>
        /// The infinity norm condition number of the matrix. The number is calculated as ||A|| * ||AInverse|| where ||.|| is the infinity norm.
        /// </summary>
        public double InfinityNormConditionNumber()
        {
            Compute();
            if (singular)
            {
                return Double.PositiveInfinity;
            }
            if (infCondition == Double.MinValue)
            {
                infCondition = matrix.lapack.ConditionNumber('I', matrix.Rows, infNorm, matrix.data, pivots);
            }

            return infCondition;
        }

        Matrix ILU.Inverse()
        {
            return Inverse();
        }

        /// <summary>
        /// Returns the inverse of the matrix.
        /// </summary>
        /// <returns>The inverse of the matrix.</returns>
        public Matrix Inverse()
        {
            DenseMatrix result = (DenseMatrix)MatrixBuilder.CreateMatrix(matrix.Rows, MatrixType.Dense);
            Inverse(result);
            return result;
        }

        /// <summary>
        /// Places the inverse of the matrix in to the result matrix.
        /// </summary>
        /// <param name="result">The matrix to place the inverse into.</param>
        public void Inverse(Matrix result)
        {
            DenseMatrix temp = new DenseMatrix(matrix);
            matrix.lapack.LUInverse(result.Rows, temp.data, pivots);
            temp.CopyTo(result);
        }

        /// <summary>
        /// Places the inverse of the matrix in to the result matrix.
        /// </summary>
        /// <param name="result">The matrix to place the inverse into.</param>
        public void Inverse(DenseMatrix result)
        {
            matrix.CopyTo(result);
            matrix.lapack.LUInverse(result.Rows, result.data, pivots);
        }

        /// <summary>
        /// Returns an <b>int</b> array indicating which rows were interchanged during factorization.
        /// Row i was interchanged with row pivots[i].</summary>
        /// <value>An <b>int</b> array indicating which rows were interchanged during factorization.</value>
        public int[] Pivots()
        {
            Compute();
            int[] ret = new int[pivots.Length];
            Buffer.BlockCopy(pivots, 0, ret, 0, pivots.Length * Constants.SizeOfInt);
            return ret;
        }

        /// <summary>
        ///  Solves a system of linear equations, <b>AX = B</b>.
        /// </summary>
        /// <param name="rightHandSide">The right hand side <see cref="Matrix"/>, <b>B</b>.</param>
        /// <returns>The left hand side <see cref="Matrix"/>, <b>X</b>.</returns>
        public Matrix Solve(Matrix rightHandSide)
        {
            Matrix result = MatrixBuilder.CreateMatrix(matrix.Rows, MatrixType.Dense);
            Solve(rightHandSide, result);
            return result;
        }

        /// <summary>
        /// Solves a system of linear equations, <b>AX = B</b>.
        /// </summary>
        /// <param name="rightHandSide">The right hand side <see cref="Matrix"/>, <b>B</b>.</param>
        /// <param name="result">The left hand side <see cref="Matrix"/>, <b>X</b>.</param>
        public void Solve(Matrix rightHandSide, Matrix result)
        {
            DenseMatrix copy = new DenseMatrix(rightHandSide);
            matrix.lapack.LUSolve(matrix.Rows, rightHandSide.Columns, matrix.data, pivots, copy.data);
            copy.CopyTo(result);
        }

        /// <summary>
        /// Solves a system of linear equations, <b>Ax = b</b>.
        /// </summary>
        /// <param name="rightHandSide">The right hand side vector, <b>b</b>.</param>
        /// <returns>The left hand side <see cref="Vector"/>, <b>x</b>.</returns>
        public Vector Solve(Vector rightHandSide)
        {
            Vector result = VectorBuilder.CreateVector(matrix.Rows, VectorType.Dense);
            Solve(rightHandSide, result);
            return result;
        }

        /// <summary>
        /// Solves a system of linear equations, <b>Ax = b</b>.
        /// </summary>
        /// <param name="rightHandSide">The right hand side <see cref="Vector"/>, <b>b</b>.</param>
        /// <param name="result">The left hand side <see cref="Vector"/>, <b>x</b>.</param>
        public void Solve(Vector rightHandSide, Vector result)
        {
            double[] temp = rightHandSide.ToArray();
            matrix.lapack.LUSolve(matrix.Rows, 1, matrix.data, pivots, temp);
            result.SetValues(temp);
        }
        #endregion Public Methods

        protected override void InternalCompute()
        {
            infNorm = matrix.InfinityNorm();
            l1Norm = matrix.L1Norm();
            pivots = new int[matrix.Rows];
            matrix.lapack.LUFactor(matrix.Rows, matrix.data, pivots);

            for (int j = 0; j < matrix.Rows; j++)
            {
                if (matrix[j, j] == 0)
                {
                    this.singular = true;
                    break;
                }
            }
        }
    }
}
